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Abstract 

We present results from numerical simulations of the Nambu— Jona-Lasinio model 
with an SU (2) <g) SU (2) chiral symmetry and iV c = 4, 8, and 16 quark colors 
at nonzero temperature. We performed the simulations by utilizing the hybrid 
Monte Carlo and hybrid Molecular Dynamics algorithms. We show that the model 
undergoes a second order phase transition. The critical exponents measured are 
consistent with the classical 3d 0(4) universality class and hence in accordance 
with the dimensional reduction scenario. We also show that the Ginzburg region is 
suppressed by a factor of 1/N C in accordance with previous analytical predictions. 
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1 Introduction 



The behavior of symmetries at finite temperature is one of the most outstanding 
problems in major areas of particle physics such as cosmology, relativistic heavy- 
ion collisions, and the quark-gluon plasma. In recent years, considerable work 
has been done on the physics of the finite temperature chiral phase transition in 
QCD and related models. Since the problem of chiral symmetry breaking and its 
restoration is intrinsically non-perturbative, the number of techniques available is 
limited. Lattice simulations have so far provided us with most of our knowledge 
about this phenomenon. Various lattice results indicate that the QCD transition is 
second order [ 1 1 . However, due to problems with finite size effects near the critical 
temperature and the proper inclusion of fermions on the lattice, these calculations 
are still not fully satisfactory. Therefore, an important role is played by model 
field theories which incorporate certain basic features of QCD and they are easier 
to deal with. Such a theory is the Nambu— Jona-Lasinio (NJL) model. This model 
was introduced in the sixties as a theory of interacting nucleons and later it was 
reformulated in terms of quark degrees of freedom. In this paper we present results 
from numerical simulations of the NJL model near the thermal phase transition. 

A compelling idea based on universality and dimensional reduction Q, is that 
in QCD with two massless quarks the physics near the transition can be described 
by the three-dimensional 0(4) -symmetric sigma model. The reasoning behind the 
dimensional reduction scenario is based on the dominance of the light degrees of 
freedom near the critical temperature T c . In the imaginary time formalism of finite 
temperature field theory, a given system is defined on a manifold S* 1 x R d , where 
d is the number of spatial dimensions. The inverse temperature is the circumfer- 
ence of S 1 . The boundary condition in the time direction is periodic for bosons 
and anti-periodic for fermions. At temperatures larger than the physical scales 
of the system, the nonzero Matsubara modes acquire a heavy mass. Therefore, 
the non-static configurations are strongly suppressed in the Boltzmann sum. If 
the original theory has fermions, all fermionic modes become massive due to the 
anti-periodicity in the temporal direction and can be integrated over. This is not to 
say that non-static modes have no effect. Rather, they generate counter-terms to 
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the three-dimensional theory. It is not yet clear whether the dimensional reduction 
is the correct description of the QCD phase transition. 

In the continuum, the NJL model is described by the Euclidean Lagrangian 
density 



where ip 0/0 are the quark (anti-quark) fields and f = (r\, r 2 , r 3 ) are the Pauli 
spin matrices, which run over the internal SU(2)j isospin symmetry. The in- 
dex % runs over N c quark colors and g 2 is the coupling constant of the four- 
fermion interaction. The model is chirally symmetric under SU(2) l <8> SU(2) R : 
ip — > {PjJJ + P R V)ip, where U and V are independent global SU (2) rotations 
and the operators Pl,r = |(1 ± 75) project onto left and right handed spinors, 
respectively. In other words the explicit flavor symmetry gives rise to Nf = 2. It 
is also invariant under U(l)v corresponding to a conserved baryon number. The 
theory becomes easier to treat, both analytically and numerically, if we introduce 
scalar and pseudo-scalar auxiliary fields denoted by a and 7?, respectively. The 
bosonised Lagrangian is 



Since the theory has no gauge fields it fails to produce the physics of confine- 
ment. The dynamic generation of fermion masses, though, brought about by the 
breaking of chiral symmetry to SU(2)j at g 2 larger than a critical g 2 is one of the 
important properties of the model. This implies the creation of both the Nambu- 
Goldstone bosons (n) and the chiral partner sigma (a), which are quark-antiquark 
bound states. The interaction strength has a mass dimension -2, implying that 
the model is non-renormalizable. The triviality of the NJL model was also de- 
mostrated numerically in J4|. This model works well in the intermediate scale 
region as an effective theory of QCD with a momentum cutoff comparable to the 
scale of chiral symmetry breaking A x sb ~ lGeV. Various authors used the NJL 
model to study properties of hadrons (for a review see 0). It was also shown that 
the model has an interesting phase diagram in the (T, ji) plane, where ji is the 




(1) 



£ = Ipi + CT + i7 5 7T • f) Ipi + ^{p 2 + 7T ■ 7?) . 



(2) 
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quark chemical potential. This phase diagram is in close agreement with predic- 
tions from other effective field theories of QCD. In the Hartree-Fock approxima- 
tion, the theory has a tricritical point in the (T, fi) plane (6| when the quark current 
mass, the coupling g, and the momentum cut-off A are chosen to have physically 
meaningful values. Recently, lattice simulations showed that the ground-state at 
high fi and low T is that of a traditional BCS superfluid Q. 

A possible loophole of the dimensional reduction scenario is that in the NJL 
model (and in QCD) the meson states are composite quark— anti-quark bound 
states whose density and size increase as T — > T c . This may imply that the 
fermionic sub-structure is apparent on physical length scales and as a consequence 
we may have a maximal violation of the bosonic character of mesons near T c . If 
this is indeed the case, then the quarks become essential degrees of freedom ir- 
respective of how heavy they are. Several authors observed |9) that in the 
large N c limit, where quarks are the only degrees of freedom, the exponents of 
the finite temperature transition are given by the Landau-Ginzburg Mean Field 
(MF) theory. The reasoning behind this result is that at finite temperature the 
zero fermionic modes are absent and therefore, contrary to a purely bosonic the- 
ory, the effect of making the temporal direction finite (1/T) is to regulate the 
infrared behavior and suppress fluctuations. However, inclusion of mesonic fluc- 
tuations in a non-perturbative scheme that takes into account 1/N C effects results 
in a first order finite temperature phase transition ifTOl . The first order transition 
discontinuity decreases with decreasing the cut-off that controls the strength of 
the mesonic fluctuations, but the nature of the transition remains unchanged even 
for very small values of the cut-off. It was also shown numerically in both the 
d = (2 + 1) Ell El and d = (3 + 1) O cases that the Z 2 -symmetric model un- 
dergoes a second order finite temperature phase transition. This transition belongs 
to the (d — 1) -dimensional Ising model universality, implying that the dimensional 
reduction scenario is a valid description. In this paper, we study the case where 
the chiral symmetry is continuum. 

In Sec. 2 we summarize the lattice formulation of the model. In Sec. 3 we 
present the observables we used in the analysis of the critical properties of the 
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theory. Finally, in Sec. 4 we present analysis of data in the vicinity of the critical 
temperature and in the low temperature broken phase. We show that the model 
undergoes a second order phase transition, which belongs to the 3d 0(4) univer- 
sality class, and that the nontrivial scaling region is suppressed by a factor of 1/N C 
in accordance with results presented in [fTTTl. 

2 Lattice Model 

The lattice action of the NJL model used in this study was the one first used in H, 

N 

S = ^2^2 (Xa(x)M[a,Tf] xy Xa(y) + ( a (x)M^[a,jr] xy ( a (y)) 

a=l xy 

2N 

+ — V(a 2 (x) + 7r(5)-vf(i)), (3) 
q z 

where x, x, £, and ( are Grassmann-valued staggered fermion fields defined on 
the sites x of a (3 + l)d Euclidean lattice. The scalar field a and the pseudoscalar 
triplet 7T are defined on the dual lattice x and the index a runs over iV staggered 
fermion species. The kinetic operator M is the usual one for four-fermion mod- 
els with staggered fermions, modified to incorporate the SU(2) ® SU(2) chiral 
symmetry and is given by 




where m is the bare fermion mass and the SU(2)i indices p, q are shown explicitly. 
The phase factors e{x) = (_i)(*o+*i+*a+*s) an d Vu = (_i)xo+...+*^-i are t h e 

remnants of the 75 and 7^ Dirac matrices on the lattice. The symbol (x, x) denotes 
the set of 16 dual sites x adjacent to the direct lattice site x. 

It can be easily shown that in the limit m — > 0, the action given by eq. © has 
a global SU(2)l ® SU(2)r 0|. In the continuum limit, fermion doubling leads 
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to a physical content for the model of 8iV fermion species - AN described by x 
and AN by £. The extra factor of 2 over the usual relation for four-dimensional 
gauge theories arises from the impossibility of even/odd partitioning in the dual- 
site approach to lattice four-fermi models. We refer to these degrees of freedom 
as colors and hence define N c = 8N. Further details about the lattice action are 
given in 0. 

We performed lattice simulations with iV c = 4, 8, and 16. The N c = 4 sim- 
ulation necessitates a fractional N = 0.5. This can be achieved using the hy- 
brid Molecular Dynamics (HMD) algorithm. More specifically, we used the "R- 
algorithm" of Gottlieb et al. lfT4l . for which the systematic error is 0(N 2 dt 2 ) 
Bl . where dt is the discrete simulation time-step. We checked that these system- 
atic errors are smaller than the statistical errors of the various observables used 
to study the critical properties of the model. On the other hand, the N c = 8 
and 16 simulations were performed by employing the exact hybrid Monte Carlo 
(HMC) algorithm. The lengths of the Hamiltonian trajectories between momen- 
tum refreshments were drawn from a Poisson distribution with mean 2.0. For the 
inversion of matrix M of eq. ©, we used the conjugate gradient method. 

The simulations were performed on asymmetric lattices with L s lattice spac- 
ings a in the three spatial dimensions, L t <C L s lattice spacings in the temporal 
direction, and volume V = L\-L t . The spatial extent must be much larger than the 
correlation length of the order parameter in order to obtain results that are close 
to the thermodynamic limit. The temperature T is given by 1/ (L t a). The temper- 
ature of the system may be altered either by varying (3 = 1/g 2 , which amounts to 
varying the lattice spacing or by varying L t . In our simulations, we kept L t fixed 
and varied the lattice spacing. 

3 Observables 

In order to study the chiral symmetry of the model we work in the chiral limit. 
Hence, we choose not to introduce a bare quark mass into the lattice action. With- 
out the benefit of this interaction, the direction of symmetry breaking changes over 
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the course of the simulation such that E = — <j(x) and IT = — J2 X K i( x ) av ~ 
erage to zero over the ensemble. It is in this way that the absence of spontaneous 
symmetry breaking on a finite lattice is enforced. Another option is to introduce an 
effective order parameter $ equal to the magnitude of the vector $ = (£, II). In 
the thermodynamic limit ($) is equal to the true order parameter (a) extrapolated 
to zero quark mass. 

The Finite Size Scaling (FSS) method 1T31 is a well-established tool for study- 
ing phase transitions. We employ this method to study the critical behavior of the 
model on lattices available to us. The correlation length £ on a finite lattice is lim- 
ited by the size of the system and consequently no true criticality can be observed. 
The dependence of a given thermodynamic observable, A, on the size L s of the 
box is singular. According to the FSS hypothesis, in the large volume limit, A is 
given by: 

A(t,L s )=L?A/»Q A (tL 1 J 1 '), (5) 

where t = {[3 C — (3)/ (3 C is the reduced temperature, v is the exponent of the cor- 
relation length, Qa is a scaling function that is not singular at zero argument, and 
Pa is the critical exponent for the quantity A. Using eq. (|5]), one can determine 
such exponents by measuring A for different values of L s . 

In the large L s limit, the FSS scaling form for the thermal average of the 
effective order parameter ($) is given by 

(*> = L^l v U{tL]/ v ). (6) 

Another quantity of interest is the susceptibility \ that is given, in the static limit 
of the fluctuation-dissipation theorem, by 

X = lim V[($ 2 ) - ($) • ($)], (7) 

where V is the lattice volume. For finite systems, the true order parameter ($) 
vanishes and for j3 > (3 C the susceptibility is given by: 

X = V(&). (8) 
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This relation should scale at criticality like Li' 1 ' . For small L s , corrections to the 
FSS relations become important. 

Furthermore, logarithmic derivatives of (&) can give additional estimates for 
v. It is easy to show that 



(&S b ) 

[ W 



(S b ) 



(9) 



where S b = 2N X](i)(° r2 (^) + ^0*0 ' ^( x )) 1S me purely bosonic part of the lattice 
action (eq. ©). Dj has a scaling relation 



D j = L 1 J»f Dj (tL] 



(10) 



4 Results 

In this section we present the analysis of numerical results in both the vicinity 
of the finite temperature phase transition and the low temperature chirally broken 
phase. The HMD R algorithm was used to simulate the NJL model with N c = 4 
whereas the HMC algorithm was used for N c = 8 and 16. 



4.1 Finite Size Scaling 

In this section we present the FSS analysis for the N c = 8 NJL model in the 
vicinity of the chiral phase transition. The lattice sizes ranged from L s = 16 to 
L s = 54 with fixed L t = 4. We used the histogram reweighting method ifTol to 
perform our study most effectively. This enabled us to calculate the observables in 
a region of couplings around the simulation coupling. We utilized the reweighting 
technique efficiently by performing simulations at slightly different couplings $ 
close to the critical coupling (3 C . Details of the simulations are listed in Tabled 
The expectation values of various observables a£ (/3) and the associated statis- 
tical errors AAi were obtained by jackknife blocking. They were then combined 
in a single expression A = A Ls {(3) to obtain, as in ifTTll . 



A 



— — + ... H - 

AAx AAi 



(AA) 2 , (11) 
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Table 1: Simulations for the FSS analysis. 





04 


Trajectories 


16 


0.5250 


220,000 


16 


0.5260 


440,000 


22 


0.5250 


200,000 


22 


0.5260 


290,000 


30 


0.5250 


100,000 


30 


0.5260 


160,000 


40 


0.5250 


65,000 


40 


0.5260 


66,000 


54 


0.5255 


40,000 


54 


0.5260 


50,000 



where AA is given by 



(AA) 2 (AAi) 2 •" (AAi) 2 ' 

In Fig. [T] we show the results of this procedure for the susceptibility x defined in 
eq. ©. 

Near (3 C the expansion of the FSS scaling relation for x, as in |[T%ll . is: 

X = a{t) + b(x)Ul' v , 
a(t) = a + + c^t 2 + ... 

b(x) = b + bxx + b 2 x 2 + ... where x = tL\ /v . (13) 

For large enough L s one can neglect a(i). We fitted the data generated on lattices 
with L s = 16, 54 to eq. (1131) by including up to the linear terms in the expan- 
sions of a(t) and b(t). The data used in the fit are shown in Fig.[T] We extracted 
C = 0.52589(9), v = 0.747(20) and 7/1/ = 1.94(5) with x 2 /DOF = 0.5. The 
values of the exponents are in good agreement with the 3d 0(4) universality ex- 
ponents v = 0.7479(90) and 7/V = 1.97(1) reported in lfT9l . The possibility 
of a first order phase transition is clearly excluded, because in such a scenario 
7/V = 3 (the number of the spatial dimensions). The high quality of this "global" 
fit is serious evidence that the reported value for the critical coupling is accurate. 
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If we neglect the a(t) terms then in order to get (3 C = 0.52572(16), v = 0.747(39) 
that are consistent with the results reported above we also have to discard the 
L s = 16, 22 data. Inclusion of the smaller lattices gives results for v and 7/1/ that 
are significantly smaller than the expected 3d 0(4) exponents. Our analysis im- 
plies that finite size effects are relatively large and therefore one has to use large 
lattices in order to reach the asymptotic FSS regime. In Fig.|2]we plot \ versus L s 
for four different values of (3. 
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Figure 1: The susceptibility x v s. ft for different lattice volumes. 

Next, we present the FSS analysis of the logarithmic derivatives Dj defined in 
eq. ©. We fitted the L s > 22 data to the linear expansion of the FSS relation of 

Dj(j3, L s ) ~ Clj l}J v + c 2j (f3 c - (3)L 2 J\ (14) 

with fixed (3 C = 0.52589 (obtained by the analysis of the susceptibility). The 
results displayed in Table |2] show that the values of v are in good agreement with 
the 3d 0(4) result v = 0.7479(90) [19|. We should also note that if we include 
the L s = 16 data in the fits, the value of v decreases and the quality of the new 
fits deteriorates as shown in Table El The outcome of the fits to D~ ~ L X J V 
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at (5 = 0.52589 is shown in Tableland plotted in Fig [21 As expected, these 
results are consistent with what we extracted from the global fits. The results of 
the analysis of the logarithmic derivatives confirm the conclusion of the previous 
paragraph that large lattices are required in order to reach the asymptotic FSS 
regime. 

In this paragraph we present the FSS analysis of the effective order parameter 
($). We fitted the L s > 22 data in the vicinity of the transition to a scaling 
function obtained from the Taylor expansion of eq. ©, up to a linear term, 

($> ~ [ Cl + c 2 (& - (3)L l J v ]L^' v . (15) 

We extracted v = 0.740(16), (3 mag /v = 0.566(8), and (3 C = 0.52593(4) with 
X 2 / DOF = 0.43. The measured value of v is in good agreement with the 3d 0(4) 
result, whereas the value of f3 mag /v deviates from the 3d 0(4) result (3 maig /v = 
0.5129(11) lfT9l . This discrepancy in the value of /5 mag might be due to higher 
order corrections to the FSS relation. In Fig. |4] we plot data used in this global 
fit versus the lattice spatial size L s at certain values of (3. After including the 
quadratic term in the Taylor expansion of the FSS relation we got results consistent 
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Table 2: Results from fits to eq. (fl4l for different ranges of L s . 



D j 


V 


X 2 /DOF 


L s = 22- 54 


D 1 


0.750(8) 


0.7 


D 2 


0.749(6) 


0.9 


D 3 


0.747(6) 


0.8 


L a = 16- 54 


D l 


0.707(4) 


1.6 


D 2 


0.714(3) 


1.7 


D 3 


0.713(3) 


1.6 



Table 3: Results from fits to Dj ~ Lj for 22 < L s < 54 at (3 = 0.52589. 





V 


XVDOF 


D 1 


0.721(23) 


0.08 


D 2 


0.721(21) 


0.15 


D 3 


0.722(20) 


0.31 
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Figure 3: D lt D 2 and D 3 vs. L s at (3 = 0.52589. 
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with those extracted from the linear expansion. Furthermore, attempts to add extra 
terms which could take into account corrections to FSS failed. The increase in the 
number of fitting parameters resulted in bad signals for the various parameters. 

4.2 Broken Phase 

We performed simulations in the broken phase in order to study the dependence of 
($) on j3 at three different values of N c = 4, 8, and 16. According to analytical and 
numerical evidence ifTTl . the non-trivial scaling region of the finite temperature 
transition in the three-dimensional Gross-Neveu model is suppressed by a factor 
of 1/ \fN~ c . The reasoning behind this result is that for finite N c , close enough to 
T c where the bosonic fluctuations become important, MF theory breaks down due 
to self-inconsistency. Hence, although the fermions do not affect the universal 
critical properties of the theory, they do influence the width of the actual critical 
region. It was also shown in IfTTl that in four-dimensions the non-trivial scaling 
region is suppressed by a factor of l/[iV c ln(L t )}. The logarithmic dependence on 
L t is an outcome of the triviality of the four-dimensional theory. 




20 25 30 35 40 45 50 55 60 



L s 

Figure 4: Effective order parameter ($) vs. L s at different values of (5. 
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Figure 5: ($) 2 vs. (3 for N c = 4, 8, and 16. The straight lines represent fits to the 
data in the mean field regions of the different iV c 's. 



0.5 



N x = 4 :-- + - 
AT = 8 x 




0.54 



Figure 6: ($) vs. /3. The continuum and dashed lines represent fits to the data in 
the 3c? 0(4) scaling regions of N c = 4, 8 repsectively. 
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We performed simulations on lattices with fixed L t = 4 and L s = 54. The 
Landau-Ginsburg MF theory predicts /3 mag = 1/2. Therefore, in the MF regions 
of different N s 's the data for ($) 2 are expected to fit well to straight lines: ($) 2 ~ 
(const. — (3) 1 . The results for ($) 2 versus j3 are presented in Fig. The iV c = 16 
data shown in Fig. |5] fit well to the MF scaling for all the values of (3. In the 
N c = 4 and 8 cases however, for sufficiently small ($), there is a crossover from 
MF scaling into the 3c? 0(4) universality. This is demonstrated more clearly in 
Fig. |6] where it is shown that the data for small values of ($), fit well to 

($> = a(/3 c (i 6) 

where (3 maig is fixed to its 3d 0(4) value of 0.3836 lfi"9l . The continuum and 
dashed lines in Fig. |6] represent fits to the data in the non-trivial scaling regions 
for N c = A and N c = 8, respectively. For N c = 4, we extracted (3 C = 0.4901(3) 
by fitting the data in the range 0.470 < (3 < 0.480. For N c = 8, we extracted 
P e = 0.5263(1) by fitting the data for couplings 0.52000 < (3 < 0.52375. In Table 
3 we show the x 2 1 DOF for different sets of data fitted to eq. (fToT l. It is clear that 
the fit quality deteriorates as we add data points at lower temperatures. Hence, we 
conclude that for N c = 4 the crossover out of the 3d 0(4) scaling region occurs 
roughly at f3 CIOSS = 0.465 with ($) croS s = 0.34 and for iV c = 8 at f3 cross = 0.5175 
with ($) croS s = 0.190. At this point we should note that in the MF region the 
order parameter is equal to the thermal fermion mass JHJ. Therefore, it acts as 
an inverse correlation length. Given that our definition of the crossover from one 
scaling region to another is a bit loose, we believe that our results are consistent 
with the analytical prediction that the scaling region is suppressed by a factor of 

i/n c ED. 

We also fitted the N c = 8 effective order parameter data for (3 =0.5200, 
0.52125, and 0.5225 to eq. dH) with the critical coupling fixed to f3 c = 0.52589 
(extracted from the FSS analysis of the susceptibility). We found /3 mag to be 
0.350(15), which is roughly two standard deviations below the 3d 0(4) result 
Anag = 0.3836(46). We believe that this small discrepancy is due to finite size 
effects which cause an increase in the value of ($). Given that (3 C is fixed, the 
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Table 4: Simulations for the FSS analysis. 



No. of points 


(3 range 


x 2 /DOF 


N c = 4 


4 


0.4725-0.4800 


0.10 


5 


0.4700-0.4800 


0.17 


6 


0.4675-0.4800 


0.58 


7 


0.4650-0.4800 


2.50 


8 


0.4625-0.4800 


3.2 


N c = 8 


3 


0.52125-0.52375 


0.08 


4 


0.52000-0.52375 


0.07 


5 


0.5175-0.52375 


7.3 



fitting fuction is forced to become more abrupt, resulting in a smaller value for 

mag- 

5 Conclusions 

We showed with relatively good precision that the NJL model with four lattice 
spacings in the temporal direction undergoes a second order phase transition. The 
critical indices were extracted by employing FSS analysis in the vicinity of T c and 
analysis of data generated on large lattices in the low temperature chirally broken 
phase. Our results are consistent with the 3d 0(4) classical spin model exponents 
and hence in accordance with the dimensional reduction scenario [2j . We summa- 
rize our results in Table |3 and compare them with those of the classical 3d 0(4) 
Heisenberg spin model |[T9ll . The values and the error bars of the exponents pre- 
sented in Table|5]take into account slightly different predictions for each exponent 
measured by different methods in this study. 

We also verified numerically a previous analytical result ifTTl that the actual 
critical region is suppressed by a factor of l/N c . Outside this region, where fluc- 
tuations are negligible, we have Landau-Ginsburg MF scaling predicted by the 
large- N c calculations. This result should be contrasted with recent numerical re- 
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Table 5: Critical exponents of the 3d 0(4) Heisenberg model obtained in [HO and 
critical exponents obtained from this study. 



3d 0(4) spin model exponents 


This study 


v 0.7479(90) 


0.75(4) 


7 1.477(18) 


1.46(6) 


/3 mag 0.3836(46) 


0.38(4) 



suits in strong coupling QCD, where the width of the actual critical region does 
not depend on iV c ll20l . 

Furthermore, it is useful to comment on the phenomenological significance 
of the parameters used in this simulation. Within the framework of the Hartree 
approximation (for details see 0), it can be easily shown that if one fixes F n to 
its experimental value of 93MeV and the constituent quark mass to a physically 
reasonable value of 400MeV then the inverse lattice spacing a -1 at the critical 
coupling [3 C = 0.5259 (for L t = 4) of the N c = 8 nonzero temperature phase 
transition is of the order of lGeV. 
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